function b_SANA=Clogit_SANA(Y_SANA,X_SANA,Z_SANA)

choice=1:13;
ns=numel(choice);

Y=sortrows(Y_SANA,[1,2]);
Y=Y(:,3:end);

X=sortrows(X_SANA,[1,2]);
X=[ones(size(X,1),1),X(:,3:end)];

Z_SANA=sortrows(Z_SANA,[1,2,3]);
Z=zeros(size(Y,1),size(Z_SANA,2),ns);
for i=4:size(Z_SANA,2)
    Z(:,i,:)=reshape(Z_SANA(:,i),ns,[])';
end
Z(:,1:3,:)=[];

%{
% X coefficients: constant bdtot nprofit teaching HHI totpop65  
b(:,2) = [-5.481; 0.00288; -0.00117; 1.676; -0.730; -0.173];
b(:,3) = [-11.33; 0.00180; 0.330; 0.659; 1.013; 0.429];
b(:,4) = [-4.771; -0.00958; 0.0883; -19.39; 0.606; 0.0512];
b(:,5) = [-2.155; -0.00994; 0.112; 1.387; -0.227; -0.226];
b(:,6) = [-10.42; 0.00370; 0.858; -0.0301; -0.936; 0.276];
b(:,7) = [-6.151; 0.000910; 1.471; 2.165; -1.783; -0.130];
b(:,8) = [-10.53; 0.00348; 0.440; 1.593; 1.322; 0.149];
b(:,9) = [-4.819; -0.0105; -0.275; -19.39; -0.556; -0.00912];
b(:,10) = [-7.879; 0.00378; 0.401; -1.257; 1.024; 0.162];
b(:,11) = [-9.054; 0.00382; 0.662; -0.818; 0.341; 0.249];
b(:,12) = [-5.693; 0.000756; 0.511; -2.667; 0.346; 0.0915];
b(:,13) = [-7.653; 0.00210; 0.545; -0.488; 0.833; 0.189];
%}

% ### UPDATED on 06/08/2020: only keep bdtot, HHI, and totpop65
%X coefficients: constant bdtot HHI totpop65  
b(:,2) = [-5.481; 0.00288; -0.730; -0.173];
b(:,3) = [-11.33; 0.00180; 1.013; 0.429];
b(:,4) = [-4.771; -0.00958;  0.606; 0.0512];
b(:,5) = [-2.155; -0.00994; -0.227; -0.226];
b(:,6) = [-10.42; 0.00370; -0.936; 0.276];
b(:,7) = [-6.151; 0.000910; -1.783; -0.130];
b(:,8) = [-10.53; 0.00348; 1.322; 0.149];
b(:,9) = [-4.819; -0.0105; -0.556; -0.00912];
b(:,10) = [-7.879; 0.00378; 1.024; 0.162];
b(:,11) = [-9.054; 0.00382; 0.341; 0.249];
b(:,12) = [-5.693; 0.000756; 0.346; 0.0915];
b(:,13) = [-7.653; 0.00210; 0.833; 0.189];



% Z coefficients: vmsh*25%Q, vmsh*25-50%Q, vmsh50-75%Q, vmsh75%Q, MKTabove25%, MKTabove50%, MKTabove75%
bz=[0.448;0.856;1.356;0.768;0.399;0.592;1.062];

bOri=[reshape(b(:,2:end),[],1);bz];
%bAns=zeros(size(bOri));

options=optimset('MaxFunEvals',200000000000,'MaxIter',1500000000000,'TolX',1e-16,'Tolfun',1e-16,'GradObj','on','DerivativeCheck','off','FinDiffType','central');
startval = bOri;
baseAlt=1;
b_SANA = fminunc('clogit',startval,options,[],Y,X,Z,baseAlt);













